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Abstract 

This work presents a method of computing Voigt functions and their derivatives, to high 
accuracy, on a uniform grid. It is based on an adaptation of Fourier-transform based convo- 
lution. The relative error of the result decreases as the fourth power of the computational ef- 
fort. Because of its use of highly vectorizable operations for its core, it can be implemented 
very efficiently in scripting language environments which provide fast vector libraries. The 
availability of the derivatives makes it suitable as a function generator for non-linear fitting 
procedures. 
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Introduction 

The computation of Voigt line profiles is an issue which has been dealt with over a 
long time in the literature yJ,l2L|3|,|4|, e.g.]. Nonetheless, it remains a computationally 
interesting problem because it is still fairly expensive to compute profiles to high 
accuracy. This paper presents a method which is fast and very simple to implement. 
It is similar to the method of [3], but capable of much higher precision for a given 
computational effort. More importantly, the method described here computes not 
only the Voigt function itself, but its derivatives with respect to both the Gaussian 
and Lorentzian widths, which are helpful for non-linear curve fitting. 

Because the method is based on Fourier transforms, and generates a grid of values 
for the function in a single operation, it is particularly suitable for (but not restricted 
to) use in scripting environments, where Fast Fourier Transforms (FFTs) and other 
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vector operations are available as part of the language library, and execute at very 
high speed. Thus, in Matlab®, Python (with Numeric or numpy), Mathematica®, 
Octave, and other similar environments, it may find particular applicability. Such 
environments provide a trade-off of speed for flexibility, and this method allows one 
to take a much smaller speed penalty than would be exacted using other algorithms 
implemented in the scripting language. 

A method such as this one, which contains no inherent approximations, has advan- 
tages over other approximate methods in that it can be adapted very easily to the 
desired level of precision of the task at hand. Selecting the density of the computa- 
tional grid and the length into the tails to which the computation is carried out sets 
the accuracy. The grid density does not affect the accuracy of the values this pro- 
duces directly, but does affect the accuracy to which interpolation may be carried 
out between the computed points. The distance into the tails to which the profile is 
computed does affect the accuracy of the step (defined below) which corrects for 
periodic boundary conditions, and it converges as the fourth power of this distance. 

As an aside on notation, most papers working with the Voigt function historically 
have defined it in terms of a single parameter, the ratio of the Lorentzian to Gaus- 
sian width. This work, presents the results in a slightly different form, which is 
more useful for direct computation. The equations below are computed in terms of 
the Lorentzian width (which I call a) and the standard deviation of the Gaussian 
distribution a. The parameter y in Dray son yj] is a/ (\/2o) in my notation. In this 
form, the transforms produce functions fully scaled and ready for direct interpola- 
tion. 



The Voigt function is a convolution of a Lorentzian profile and a Gaussian, 



and can be easily written down in Fourier transform space using the convolution 
theorem: 



Also, of great importance to using this in fitting procedures, the derivatives of this 
function with respect to its parameters can be computed: 



Theory 




(1) 



V(oc,g; k) =exp [-o 2 k 2 /2-a\k\] 



(2) 
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and 

ay 2 ~ 

— = -ok 2 V (4) 
da 

and since the differentiation in these cases commutes with the transform, these are 
the transforms of the appropriate derivatives of the function itself. 

Note that, since this transform method generates functions with a fixed area, these 
are the derivatives with respect to the widths at fixed area, rather than at fixed 
amplitude. This implies that fitting carried out using functions computed this way 
is most appropriately done using a, a, and area as parameters. 

This result is exactly correct in the full Fourier transform space. For practical com- 
putation, though, one wishes to reduce this into something which is computed 
rapidly on a discrete lattice using Fast Fourier Transform (FFT) techniques, and 
then interpolated between the lattice points. The difference between the full contin- 
uous transform and the discrete transform is, of course, that the function produced 
by a discrete transform is periodic. In effect, by discretely sampling the series in 
Fourier space, one is computing, instead of the exact convolution, a closely re- 
lated function which has had periodic boundary conditions applied. This affects the 
shape of the tails of the distribution, but in a way which is fairly easily fixed. 

First, when doing the discrete transforms, it is necessary to decide how far out in 
fc-space it is necessary to have data. In general, one wants to assure the function 
is nicely band-limited, which means no significant power exists at the highest k. 
Practically speaking, setting the argument of the exponential in eq.[2]to something 
like -25 at the boundary means the highest frequency component is e~ 25 or about 
10 11 of the DC component. To achieve this, define the absolute value of the log 
of the tolerance to be y (25 for the example here) and solve o 2 k 2 /2 + <x\k\ = y. 
This is a simple quadratic equation, but because one doesn't know in advance how 
dominant the relative terms are, one should solve it with a bit of care. The stable 
quadratic solution presented in Numerical Recipes JH] can be adapted to be 

r (5) 



■~mux I — ^ =■ 

a + -y/a 2 + 2ya 2 

This is simpler than the full solution in [5] since the signs of both a and o are 
known in advance. 

Now, note that the periodic solution is really an infinite comb of functions shaped 
like the desired one, added together. Since, beyond a few a from the center, the 
function is close to Lorentzian, one has really computed the desired function plus 
an infinite series of offset simple Lorentzians: 

V (a, a; x) = V act (a, a; x) + £ £ (jc _J )2 + a2 (6) 
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where A is the period of the function. However, the infinite sum can be computed 
analytically. It is: 

cc^ 1 sinh 2 ^ /a\ 1 



- = a im t (7) 

n ^ Q (x-nA) 2 + a 2 A(cosh^-cos^) \nJ x 2 + a 2 

de _ (a 2 -x 2 ) 27icosh(^) 27tsinh 2 (^) 

3a " n (a 2 +x 2 ) 2 + A 2 (cosh ^ - cos ?f) ~ A 2 ( cosh 2™ _ cos 2^ 2 (8) 



Since the derivative with respect to a is very localized, and falls to zero rapidly at 
the boundaries, no correction is needed for it. 

Although these equations look computationally intensive, they are not so at all. 
Note that the cosh and sinh terms are of a constant, and not evaluated at each point. 
Also, for the usual Fourier-space case, A = 2x max so the cos term is just an evalu- 
ation of the cosine from —n to % on the same grid the rest of the function will be 
evaluated. If the Voigt function is to be evaluated for many different a, a pairs (as 
is the case in fitting routines), but always on a grid with a fixed number of points, 
this cosine only gets evaluated once, too, and can be cached for reuse. Also, the 
correction and its derivative with respect to a share most of their terms in common, 
so this correction is really a simple algebraic adjustment to the raw function table. 

The correction term in eq. [7J is an approximation based on all the other nearby 
peaks being entirely Lorentzian, and works well. However, it can be improved by 
a scaling argument, which works significantly better. The non-Lorentzian nature 
of the correction is due to the convolution of the Gaussian with the curvature of 
the Lorentzian causing a slight widening even on the tails. Note that convolution 
of a function with a Gaussian only affects even derivatives (by symmetry), so the 
second derivative term is the lowest order this could affect. Also note that this effect 
is getting bigger as one approaches the next peak over (the edge of the boundary). 
Thus, one can try a correction of multiplying the right hand side of eq. [TJby 1 + ax 2 
where a is to be determined. Much of the structure of a can be obtained by scaling, 
and it should be a oc a 2 /A 4 . Empirical testing has shown that a constant of 32 
appears optimal, so an improvement on eq. [TJis: 

_ I" sinh 2 ^ /cc\ 1 

fc ~ A^osh^-cos 2 ^) ~ \ n) x 2 + a 2 



32oV 
A 4- 



(9) 



This improves the original correction by almost an order of magnitude in the peak 
error at the bounds of the interval, and the RMS error is reduced by about a factor 
of 5 for most test cases. 

Note that this expression also provides an error estimate for the calculation, which 
can be used to determine an appropriate value for A. Assuming the error is of the 
same order as the final correction term, which should be conservative, one can 
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proceed to evaluate it at the boundary of the interval, where x = A/2, and evaluating 
eq. |7] assuming A ^> a, the error estimate is then 



relative err ~ (small coefficient) x 2o 2 /A 2 



(10) 



Note that, although this is 1/A 2 , this is the relative error. The function itself 
is decreasing at the boundary as 1/A 2 so the absolute error scales with 1/A 4 , as 
expected. 



The road not taken 

There is another way one could consider carrying out this computation, which looks 
elegant and easy from the outset, but actually is computationally much more expen- 
sive. I will outline it here as a warning to others. 

Instead of fixing the periodicity error by adding on the correction of eq. [7] or eq. 
[9l one might be tempted to fix the problem in advance, before the transforms are 
carried out from fc-space to real space. The obvious solution is to try to compute 
the transform of the difference between the Voigt function and a pure Lorentzian, 
and then add the pure Lorentzian back in afterward, not as an infinite sum as in the 
correction equations, but just as a single copy. One would compute 



and then transform this, and add back on the Lorentzian which was subtracted. 
This turns out to be computationally very inefficient, though. When computing the 
transforms, one wants a cleanly band-limited function in £-space, with the power in 
the highest frequency channels vanishing rapidly. In the case of eq.[2l this is clearly 
the case, since the —o 2 k 2 /2 term makes the exponential disappear relatively rapidly 
even for fairly modest values of o and k. In the case of eq. [EH though, % only 
vanishes as fast as the exp[— a\k\] term, which falls off much more slowly. Thus, 
one has to carry out the transforms to much higher values of k to get convergence. 
This turns out practically to be a huge penalty. Even in the case of o ~ k, it requires 
about a few times more terms, and in the case o k, it is much worse, since the 
extremely rapid falloff of the Gaussian in &-space allows one to sample only quite 
small values of k to get very good performance. 



Application 

The most probable way the author sees these gridded functions being used is to load 
cubic spline interpolation tables to generate values at points which may not lie on 



Vb(a,a; k) = (exp [-a 2 fc 2 /2] - l) exp[-a|&|] 



(11) 
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the grid. This way, one can compute the Fourier transforms on grids of sizes con- 
venient for FFT algorithms (often, powers of two, but using, e.g. FFTW [6], many 
other grid sizes can be conveniently transformed), and then use the interpolator to 
fill in the values desired. Because of both the shape of these functions and the wide 
dynamic range they typically encompass, it is likely that it will be useful to interpo- 
late the logarithm of the function. Especially if a 3> a, so the center looks Gaussian, 
log interpolation is extremely beneficial, since the logarithm of the Gaussian part is 
just parabolic, and exactly interpolable by a cubic spline interpolator. 

In fitting work, one more derivative is needed than the ones computed in eqs. |3]and 
IH which is the derivative with respect to the coordinate itself, which is needed to 
solve for the center of a peak. Although this could easily be computed in Fourier 
transform space by inverse transforming ikV(k), it also falls directly out of the 
cubic spline interpolation. The value of a splined function at a point is computed 
(see the discussion of splining in [5Q for notation) as 

h = x\ — xq 

X\— X 

a= — 

b= \ — a 

y = ayo + b yi + j [(a 3 - a) yg + (b 3 - b) y'{] (12) 

then the first derivative is 

y = ^ + ^[(3^ 2 -lM-(3a 2 -l)ja (13) 

This is the method preferred by the author of this work for this derivative. 



Error Analysis 



Figure Q] shows sample functions computed by this method, and the relative er- 
rors associated with this computation. These were computed using the extra error 
correction of eq. |9l The scaling of the errors is fairly easy to compute from the 
underlying equations. In general, the errors scale as A~ 4 when a and a are held 
constant. For most practical applications, it is likely that the need to compute the 
tails far enough from the center that the entire spectrum is covered by the calcula- 
tion results in the tails being calculated sufficiently far out that the accuracy is not a 
concern. In Figure [H the curve (5) shows the result when A = 80a (tails computed 
to 40a), and even in this case the peak relative error is 10~ 4 (in a part of the curve 
off the graph). Most practical cases are likely to need the tails much farther out than 
this, resulting in the accuracy automatically being better than this. As an example, 
though, of pushing the computation to very narrow tails, the curve (4) shows the 
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results for the tails being computed to only 10a. Even in this case, the relative error 
only exceeds 10 -3 at the very edge of the domain. 

Note that these are computed with a = 1 and a being varied. Since the shape of 
the Voigt function really only depends on a/a (related to the usual y parameter), 
this completely specifies the shape of the function except for an overall width scale. 
Also, note that the cases shown are all for a > a; the case of small a is a limiting 
case, and the errors look much as they do for a = a. 

Conclusion 

Because of relatively slow convergence, simple FFT-based convolution has not 
fared well in the Voigt-function computation arena. Nonetheless, this method has 
always had an advantage in simplicity and vectorizability. Also, it is trivial to get 
the derivatives of the Voigt function with respect to its width parameters from 
transform-based methods. This makes this algorithm most useful for cases in which 
line strengths, widths, and positions are all variable. The convergence enhancement 
can also be easily extended to line shapes in which a non-Gaussian function is con- 
volved with a Lorentzian. 

By combining the traditional transform-based method with a convergence-enhanc- 
ing operation, the result is a method which is fast, accurate, and extremely easy 
to implement. It should find particular application for fitting work carried out in 
many widely used scripting languages, in which fast vector operations often make 
computation of tables of function values an efficient process. As an example, on 
my 1 GHz laptop computer, it takes 8 milliseconds to compute a 2048 point grid 
of the function and its two derivatives using this method, in the Python language. 
Even in compiled languages, though, this may be highly adaptable to fast work on 
any operating system and machine which provides good vector operation and FFT 
support. 

A more detailed speed comparison of this to other algorithms, in other languages, 
with differing numbers of points, and differing accuracy requirements, is difficult, 
since every one of these parameters affects the speed of one algorithm relative to 
another. Compared to explicit, point by point, computation of the Voigt function 
in a scripting language (which was done, for Fig. [TJ by direct convolution) it can 
be 2 orders of magnitude faster. Compared to the unenhanced, transform-based 
algorithm of Karp [3D, this provides much more accuracy for the same speed, or 
many times better speed for the same accuracy, where 'many' can often be a factor 
of 10. 
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Fig. 1. Plot of Voigt Functions with a = 1 evaluated by this method, and the relative errors 
(as computed by adaptive numerical integration of the convolution at each point). 
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